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ABSTRACT 

This paper presents a number of variations on the Davenport algorithm for in-flight gyroscope recalibration, or 
first-order initial calibration, specifically tailored for use with a minimum amount of satellite telemetry data. 
Central to one of the techniques described is the use of onboard integration of gyroscope data together with a 
detailed model of scheduled satellite slew profiles. Methods are presented for determining adjustments to either 
parameters for the standard linear model (i.e., a drift rate bias vector and/or a scale factor / alignment 
transformation matrix) or individual gyroscope scale parameters, both linear and nonlinear, in cases where the 
alignments are well known. The results of applying the methods in an analysis of the temporal evolution and 
nonlinear response of the gyroscopes installed on the Hubble Space Telescope following its first servicing mission 
are discussed. The two effects, when working coherently, have been found to result in slew errors of almost 
1 arcsecond per degree. Procedures for selecting optimal operational gyroscope parameters subject to the 
constraint of using a linear model are discussed. 

Introduction 

Reference 1 presents a derivation of the Davenport 
gyroscope calibration algorithm, which has been 
used for the in-flight calibration of gyroscopes for a 
number of spacecraft missions, including those of 
the High Energy Astrophysics Observatories and the 
Hubble Space Telescope (HST). As usually 
implemented, and, in particular, as implemented for 
the HST mission (Reference 2), the algorithm 
assumes that the user has available for use in the 
calibration process a continuous and complete set of 
gyroscope data extending from an initial to a final 
spacecraft attitude (as determined by independent 
reference sensors) for an adequately large number of 
maneuvers. Empirically, we find that this constraint 
causes gyroscope scale factor and alignment 
calibration to be one of the more labor- and data- 
intensive activities needed in support of mission 
operations. Fortunately, we also have found that the 
scale factor and alignment parameters for the 
gyroscopes used for the HST mission are fairly 
stable; calibration is usually required only following 
initial deployment of gyroscopes (i.e., following 
HST’s initial deployment in April 1990, activation 
of reserve gyroscopes in response to gyroscope 
failures, and installation of new gyroscopes during 
the first HST servicing mission in December 1993). 


Although the HST gyroscopes are “fairly stable,” a 
performance analysis conducted in September 1995 
(Reference 3) has indicated that in the 18 months 
following the first HST servicing mission, the 
gyroscope response has changed systematically, the 
errors being most manifest in negative yaw 
maneuvers wherein systematic errors of roughly 
0.8 arcsecond per degree occur. 

Given this recent experience with the HST 
gyroscopes, we have found it desirable to develop 
an algorithm that permits recalibration of the 
gyroscopes, at least to first order in the change 
parameters, using a data set that is both much 
reduced in volume and readily available during 
normal mission operations. We also have found it 
useful to extend the algorithm to allow study of both 
isolated and nonlinear scale corrections. The 
algorithm that we present here requires as input 
from telemetry only the attitude error measurements 
determined by the onboard computer (OBC) 
pointing control subsystem following large vehicle 
maneuvers. All other required input can be obtained 
from the schedule of commanded maneuvers and the 
spacecraft parameters characterizing those 
maneuvers. 
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The body of this paper is divided into six sections, 
excluding this introduction. These include (1) back- 
ground on the basics of the Davenport algorithm, 
(2) a reformulation taking advantage of OBC inte- 
gration of gyroscope data and modeling of planned 
maneuver profiles, (3) some comments on the cali- 
bration of gyro bias, (4) an extension to both 
isolated and nonlinear scale factor corrections, (5) a 
discussi on of selection of measurement weights to 
be used in the algorithms, and (6) an application of 
the algorithm to data accumulated for the gyro- 
scopes used for the HST mission. 

The Davenport gyroscope calibration algorithm, as 
well as the variations of it discussed in this paper, 
are envisioned as applied in a batch mode least- 
squares algorithm. Batch mode processing is strictly 
appropriate only if the time scale for collection of 
the calibration data is short compared with the time 
scale for any variation that may apply to the state 
vector parameters. Empirically, in the case of the 
gyroscopes used on HST, we have found the scale 
factor and alignment parameters sufficiently stable 
that a batch mode approach for their calibration is 
operationally viable. In cases where this fails to be 
true, reformulating the calibration equations 
presented here in terms of a Kalman filter (e.g., 
Reference 4) should be considered. 

Section 1 - Background on the Basics of the 
Davenport Algorithm 

Reference 1 presents the gyroscope calibration 
algorithm that is used in the HST mission for the 
calibration of scale factors, alignments, and biases 
of the gyroscopes when one or more gyroscopes are 
first activated for operational use. The basic 
equations are as follows. Consider a satellite gyro 
system consisting of N 0 single-axis gyroscopes. In 
response to some angular motion of the satellite, the 
output response column matrix of gyro counts, C, 
consists of the N 0 individual gyro readings. The 
response vector is translated into a measured angular 
velocity, Q M , in the spacecraft frame via 

q m = g () c-d 0 (1) 

where G 0 is the 3xN 0 scale factor / alignment matrix, 
and D„ is the gyro system drift rate bias expressed in 
the spacecraft frame. The goal of the algorithm is to 
determine correction matrices m and d that may be 
applied to G„ and D 0 so that a modified equation (1) 
will yield the true angular rate, £2, as indicated in 
equations (2a) - (2c). 


G = (/, + m) G 0 

D = (/, + m) D # + d (2b) 

a = G C - D = (/, + m) Q„ - d (2c) 

where /, is the 3x3 identity matrix. Gyroscope 
miscalibration information is sampled through any 
given maneuver via the error quaternion 

8Q - Q,Q 0 ' ( 3 ) 

where Q R represents the true vehicle rotation as 
determined from reference star measurements, and 
Q 0 represents the vehicle rotation inferred from the 
gyroscope measurements. As discussed in 
Reference 1, 8Q represents a rotation from the 
gyro-inferred to the true postmaneuver attitude, 
expressed in the premaneuver reference frame. The 
information content of 6Q is related to m and d via 
the sensitivity equation 

Z = l/2jr(Q-Q M )dt (4a) 

= 1/2 j T ( m Q M - d ) dt (4b) 

where Z is the vector component of 5Q, T is the 
matrix that transforms vectors to premaneuver 
spacecraft coordinates, and the time integral is over 
the whole maneuver. Because equation (4b) is 
linear in m and d, it can be used as the basis for a 
linear least-squares algorithm to provide estimates 
for m and d. If a solution for all 12 correction terms 
is needed, at least 4 independent “maneuvers” are 
required to perform the calibration. The maneuvers 
must provide a reasonable sample of pitch, roll, and 
yaw variation, as well as an independent sample for 
bias determination; the latter is permitted to be a 
period of essentially constant attitude. 

Although the information content is unchanged, it is 
often more convenient to reexpress Z in terms of an 
error vector, £, representing the rotation from the 
true postmaneuver attitude to the intended (and 
gyro-inferred, assuming closed-loop control) post- 
maneuver attitude, i.e., the rotation that the 
spacecraft must perform after it determines its post- 
maneuver error. The vector £ is related to Z and 
(in, d) via 

C = - IV 1 z 

= - l/2r T ‘Jr(m£2 M -d)dt (5) 
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where i represents the maneuver duration time, and 
its use as a subscript on T x means that T is to be 
evaluated at the maneuver end-time. The matrix 7Y 1 
(which equals T x ) is thus the premaneuver to post- 
maneuver reference frame transformation matrix. 

Section 2 - Use of OBC Gyro Data Integration 
and Model Maneuver Profiles 

As discussed in Reference 1 , equation (5) is accurate 
only to first order in m and d, implying that the 
associated least-squares algorithm is intrinsically 
iterative. The matrix terms C. T, and £i M must be 
reevaluated on each iteration. Multiple iterations 
can only be applied if a complete set of gyroscope 
data from throughout each of the calibration 
maneuvers is available. In this section we discuss a 
procedure that excludes the possibility of multiple 
iterations, the gain being a drastic reduction in the 
total volume of data required to perform the calibra- 
tion. This can be significant if either (1) the sheer 
volume of data for frequent, normal calibrations 
becomes unwieldy or (2) the standard telemetry 
format used does not contain an adequately dense 
sampling of gyro data for accurate integration. 

If calibration needs are adequately met via a first 
order correction, it is possible to implement an 
algorithm with drastically lower data requirements. 
Integration of the full set of gyroscope data is 
required at two points in the use of equation (5): 
first in the determination of Q 0 for the construction 
of £, and then in the time integral over (Tmil u ). 
Ground processing of gyro data to determine Q 0 can 
be eliminated if the spacecraft OBC maintains and 
transmits an estimate of the spacecraft attitude based 
solely on gyroscope data, at least through the time 
period between the accumulation of star sensor riata 
for pre- and postmaneuver definitive a tt itude 
estimation. Sampling the OBC’s pre- and postma- 
neuvcr attitude estimates then allows construction of 
Q 0 as the connecting eigenvector rotation between 
the two. Ground processing of the gyro data for use 
in the integral over (TTnQJ can be eliminated if a 
sufficiently precise model of the maneuver profile is 
available. This follows because, to first order in the 
correction terms, equation (5) is unchanged if mil u 
is replaced with mil p , il p being the planned angular 
velocity as a function of time based on spacecraft 
design parameters. We make the latter substitution 
in what follows. 

The simplifications noted in the preceding para- 
graph allow the elimination of all ground processing 
of the raw gyro data. The elimination of ground 


processing of the reference star data may also be 
possible, although this results in a smaller gain. For 
many satellites, the OBC generates an attitude error 
estimate based upon postmaneuver reference star 
measurements and uses this estimate to generate an 
error nulling maneuver. If the vehicle attitude is 
maintained accurately by the onboard pointing 
control system during the periods between maneu- 
vers, the postmaneuver error nulling maneuver will 
correspond to the error vector £ required for our 
analysis. If this error vector is included in telem- 
etry, no other spacecraft data are required. 

We assume finally that each maneuver is a pure 
eigenaxis maneuver. This allows the analysis to be 
done in a coordinate system, designated here with a 
prime Q, in which the x'-axis is aligned along the 
maneuver axis. Expressed in the primed frame, 
equation (5) becomes 

* (2 C) =-r x T Jr(m'fl' F -d')dt (6a) 

10 O' 

T' = 0 cos[0(t)] -sin[0(t)] 

0 sin[0(t)] cos[0(t)] 


= RTR=RTR t (6b) 

m - R m R ] = R m R 1 (6c) 

d' = #fd (6d) 

Q' P = [ <n(t), 0, 0] T = RQ, (6e) 


where R is the transformation matrix that converts 
premaneuver spacecraft coordinates to the prema- 
neuver primed frame, 0(t) is the maneuver angle as 
a function of time, and to = d0/dt. The form of 0(t) 
will depend upon the total maneuver angle, 0, and 
design parameters governing the execution of 
maneuvers. To first order, R may be based on the 
planned maneuver quaternion, Q p . The eigenvector 
and rotation angle set, {r, tp), defining the quater- 
nion representation of R is constructed from the 
spacecraft frame Q, eigenvector, q, and the space- 
craft frame standard unit vectors, (x, y, z), using 

r = (xxq)/lxxql 

= ( -yq, + zq 2 ) / (q 2 J + q, J ) l/1 (7a) 

<p = cos '( x - q ) = cos ‘( q, ) (7b) 



The simple forms of equations (6b) and (6e) allow 


equation (6a) to be reexpressed as 

(2 C) = -T\ T [K l [m'] l - K„<n (8a) 
K k 0 O' 

K k m 0 K a -K* (8b) 

.0 K* 

= j (0 k dt (8c) 

K* = J cos(9) of dt (8d) 

K* = j sin(0) co k dt (8e) 


where [m\ indicates the column matrix formed 
from the first column of m‘. Note that the elements 
of Jf, are analytic, i.e., K, = ©, K c , = sin(0), and 
K„ = [l-cos(0)l, whereas K 0 is equal to the 
maneuver duration, t. The functional form of 0(t) 
enters only via K a and K l0 . 

The multiplication of T\ J into K t and K„ in equa- 
tion (8a), together with an application of the sine 
and cosine laws for two angle sums, produces 


R (2 0 = -[*■*, - **„<!'] (9a) 

K k 0 O' 

K* k = 0 K% (9b) 

0 -K** 

K*^ a j cos(0-0) of dt (9c) 

K*. = j sin(0-0) of dt (9d) 


Because A, depends only on 0 and not the form of 
0(t), it can be shown that K*, = A", 1 . This 
relationship holds for K* 0 as well (actually, for all k) 
if to(t) is an even function of time about the 
maneuver midpoint. This constraint, which is fairly 
standard for spacecraft maneuver profiles, also 
yields the following convenient equations for 
and K rt (expressed for general k): 

(10a) 
(10b) 


We now need to transform equation (9a) back into 
the spacecraft frame so as to have £ related to m and 
d rather than to m’ and d'. Defining m as the 9-by-l 
column matrix [[m], T , [m] 2 t , [m]j T f and using equa- 
tions (6b) and (6c), we can rewrite equation (9a) as 

2 £ = - R r K*, Bm + R t K* a R d (11a) 

b.,. 30.,, - R^R,. (Hb) 

where equation (11b) defines the elements of the 
3-by-9 matrix B. Equation (11a) is our new least- 
squares algorithm sensitivity equation. Its use 
removes the need for an integration of the gyro 
telemetry data. The only required time integrations 
are for K*^ and K* rt , or, more simply, F o (0) if the 
symmetry constraint on to(t) is applied. Appexdix A 
presents a specific, fairly common maneuver profile 
usable in the latter evaluation. 

Section 3 - Bias-Only Calibration Assuming Fixed 
Scale and Alignment 

We consider now the application of the algorithm of 
Section 2 to a bias-only calibration. This begins 
with the constraining assumption m - 0. This 
constraint is reasonable for many operational 
scenarios; empirically, it has been found that 
spacecraft gyro biases can change significantly 
within as little as a single day, whereas time scales 
for scale factor and alignment are considerably 
longer. For this situation, equation (11a) reduces to 

2C = R J K* 0 Rd (12) 

Two data gathering scenarios are of possible interest 
for this calibration. For the first scenario, the 
spacecraft is held at constant, or nearly constant, 
attitude over the time period of interest. “Nearly 
constant” in this context means that the magnitude 
of any net maneuver angle must be smaller than the 
product 8dAt, where 8d is the maximum permitted 
error in the estimate for d, and At is the time period 
between two reference attitude measurements. For 
this case, [R T K* a R] reduces to /,At, and equa- 
tion (12) becomes 

d = 2 £ / At (13) 

We have used At rather than i here because there is 
no scheduled or executed maneuver for which we 
can evaluate i(0). The vector £ may be constructed 
from separate initial and final reference star 
measurements or from an OBC -determined attitude 


K* = cos(0/2) F k (0) 

K* = sin(0/2) F k (0) 

F k (0) = | cos [0(t) - 0/2] of dt 


(10c) 
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error at the end of the time period if the spacecraft 
applied an attitude correction at the start of the 
period. In the latter case, care must be taken to 
ensure that the onboard attitude propagation across 
the time period involved the use of gyroscope data 
only, i.e., no control-law feedback based on 
reference star data. 

The second data gathering scenario uses the proce- 
dures outlined in Section 2 applied to equation (12). 
The potential operational advantage of using this 
approach arises if a set of dual mode gyroscopes, 
i.e., sensors with high-rate and low-rate modes, is 
being used — the high-rate mode being used during 
large maneuvers to allow greater dynamic range, 
and the low-rate mode during periods of near- 
constant attitude to allow greater precision. For 
such gyroscopes, equation (13) can be used to 
calibrate the high-rate mode bias only if the gyro- 
scopes are commanded to remain in high-rate mode 
during the calibration period, implying that 
dedicated spacecraft time would be required for the 
calibration. In contrast, the use of equation (12) 
would allow relatively frequent high-rate mode bias 
calibrations based on serendipitous maneuvers. 

A caveat pertains here - one of relevance to the 
next section. An estimate of the bias based on equa- 
tion (12) applied to a single maneuver may fail to be 
good if the estimate for the scale factor / alignment 
matrix G is insufficiently accurate, because of either 
poor initial calibration or an actual change in the 
gyroscope parameters since the time of calibration. 
The effect of errors in estimates for linear scale 
factors would tend to cancel each other in the 
estimate for d if the least-squares fit is performed 
using an ensemble of randomly directed maneuvers 
or paired sets of oppositely directed maneuvers. 
Taking advantage of this fact to reduce the influence 
of possible scale factor errors may be desirable. 

Section 4 - Isolated and Nonlinear Scale Factor 
Calibration 

The original Davenport algorithm combines the 
observable aspects of alignment and scale factor 
changes into the single change matrix m. It also 
assumes that gyroscope response is purely linear. 
We have found it useful to be able to study the 
indiviual gyroscope response curves, with respect to 
both nonlinear corrections as well as temporal 
variations of the dominant (i.e., linear) terms. In 
this section we discuss an extension of the algorithm 
presented in Section 2 designed for this purpose. 


To keep the initial discussion simple, we will 
assume that the state vector for our problem is 
restricted to scale factor adjustments. Specifically, 
we assume that the gyro alignments are well known 
and fixed, and that the adjustment to the operational 
drift rate bias is restricted to that associated with 
scale factor corrections (i.e., the “wiD 0 ” term in 
equation (2b)), with no intrinsic bias changes 
permitted (i.e., no d term in the state vector). We 
will relax both of these simplifications eventually. 
We will, however, not attempt to model ongoing 
temporal changes in the drift rate bias that occur 
during the time period over which the calibration 
data are accumulated. Given that bias changes 
occur relatively rapidly and are likely to be signif- 
icant over the data accumulation time period, this 
last simplification may at first glance seem inappro- 
priate. If, however, a good estimate of the changing 
bias vector associated with the operational align- 
ment and scale calibration is maintained throughout 
the period of data accumulation (using the methods 
of Section 3), the effect of the bias will have been 
removed on an ongoing basis. From this perspec- 
tive, we see that we are not actually neglecting the 
changing bias; rather, the bias effects have been 
precorrected as part of ongoing operations. 

To allow for nonlinear scale effects, we assume a 
model for gyroscope responsivity of the form 

C Q = S n [Q^ + I Snkgk (Q m ,)] (14) 

where subscript n indicates the nth gyroscope, C D is 
the resultant gyro reading, S D is the nominal (or 
current best estimated) gyro scale factor, ft m4i is the 
spacecraft angular rate projected onto the gyro input 
axis, and the summation over k represents a set of 
small corrections to the predominantly linear 
relationship between ft^ and C D . The parameters s^ 
are correction coefficients applied to the functions 
The latter can be any convenient set of 
functions, subject only to the constraint that the 
same set of functions be used for all of the 
gyroscopes. To minimize the eventual size of the 
least-squares state vector, the functions should be 
selected so that a good fit can be found with as few 
correction functions as possible. For our HST 
analysis, we have found it convenient to use two: 
g,(ft) = ft and g 2 (ft) = g,(ft) = Iftl. In this model, 
s nl represents an average linear correction, and s n2 
represents the difference between scale factors for 
positive and negative maneuvers. 
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We assume next that an acceptably accurate inverse 
to equation (14) can be written in the form 

= ca + I°-&(CA) 05) 

In principle, each cr^ is a function of the full set 
...}. However, if the sum Is^Q^) and 
all of its individual terms are small relative to £2^, 
and if the correction functions g k vary continuously, 
then a* = -s^ to first order in the correction terms 
for all n and k values. We will be using this 
approximation in what follows. 

For notational compactness, equation (14) can be 
rewritten as 

C = S(A£2 + £s k [*(AQ)]) (16) 

where C is a N 0 -by-l column matrix (the C/s), S 
and s t are N 0 -by-N 0 diagonal matrices (the S 0 *s and 
s^’s), A is the N 0 -by-3 matrix of gyro input axis unit 
direction vectors, and the symbol [**'] is defined 
such that 

p- V] « [fcOU, &(V 2 ), g k (V N )l T (17) 

for any N-by-1 column matrix V. We also need a 
matrix version of equation (15) that gives Q M as a 
function of C. If N 0 exceeds 3, our equation must 
include a weighting scheme for how the gyro data 
are to be combined in forming £2,*. We use the 
following equation: 

4* = [A T A] , A T Q ta 

= K (S‘C + Ic k r (S‘C)]) (18) 

where G„ is the N 0 -by-l matrix formed from the 
various G*, estimates, X = [A t A]'A t , and o k is an 
N 0 -by-N 0 diagonal matrix (the a to ’s). By using 
equation (18) as the mechanism for constructing G M 
from C, we have selected a convention whereby 
equal weight is given to each of the components of 
G to . This is a change from the more typical 
convention in constructing the matrix G 0 for 
equation (1) whereby equal weight is given to each 
component of C. 

At this point we should clarify notation a bit in 
preparation for constructing the least-squares 
algorithm for a recalibration of the s* coefficients. 
Equation (16) should be viewed as applying the 
true s* values; it represents the actual response of 


the sensors. In contrast, equation (18) represents the 
users interpretation of the counts; thus the a* are 
functions of (s„ 10 , s^, ...), where the subscript 0 
indicates current estimate. The “small correction 
terms” approximation thus leads to a* = -s* 0 . The 
least-squares state vector will be the set { 8s* } for 
all n and k, where 8s* = s* - s* 0 . 

To proceed with a formulation of an extended least- 
squares algorithm based on equation (4a), we 
require an expression for ( G - Q M ) linear in the 
correction terms 8s k . Combining equations (16) and 
(18) yields 

G„ = X {AO + Zs k [“ AO] 

+ Z a k r (AO + Z s K r AG])] ) (19) 

The assumptions that the g k functions vary 
continuously and that all of the s k and <r k elements 
are small imply that terms of the form 
c k [“ (AO + Zr, [** AO])] are equal to c k [“ AO] to 
first order. Using this simplification and setting 
0 k = - s w yields 

(0-0*) = - X X { 8s k [“ ' AO] } (20) 

To be able to follow our analytic maneuver model 
approach as developed in Section 2, we insert equa- 
tion (20) into equation (4a) and apply appropriate 
transformations to the “primed” reference frame. 
The resulting sensitivity equation is 

2C = Jf T r l T Jr(«xs{8s k [“[A/f T ]oi])dt 

= jr T X{r, T fr (* x r • [wA* T ] 1 ] D )dt}[& k ] c (2D 

where [“ [<ft4Jl T ],] D indicates a diagonal matrix 
formed from the elements of [“ ' [aiA/f 7 ],], and [8s k ] c 
indicates the column matrix formed from the 
diagonal elements of 8s k (recall that &s k is a diagonal 
matrix). If we impose the additional constraint on 
each g* that it satisfy the commutivity relation 
&(ab) = g k (a)g k (b), equation (21) can be written in 
the convenient form 

2 C = R 7 1K* 0 (R X[“ [8s k ] c (22) 

where the If** matrices are defined analogously to 
the K* t matrices discussed in Section 2, with g,(<o) 
replacing af in defining the required components. 
For the case of {g,(0)sO; g 2 (Q) = g,,(0) = IGI), 
equation (22) becomes 
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2 C = « <[[A* T ],] D [8sJ c 

+ t M [AJf T ],] D [&J C ) (23) 

where the symbol [" ] in the last term implies that 
the absolute value operation is applied to all of the 
elements of [A/? 7 ],. The matrix applies to the 
last term with no adjustments because co is by 
definition positive in the primed reference frame. 

For each maneuver used in the calibration process, 
equation (22) provides three linear equations in the 
N o ,k m« unknowns {Ss^}. To get proper visibility for 
accurately measuring all of the (Ss^J elements, a 
range of both positive and negative maneuvers in all 
of the pitch, roll, and yaw directions must be 
sampled. With an appropriately large number of 
maneuvers sampled, equation (22) can be used as 
the basis for a standard least-squares algorithm to 
determine estimates for the correction terms. 

As with the original Davenport approach to the 
calibration problem, adjustments to the scale factor 
calibration imply an associated adjustment to the 
current estimate for the drift rate bias vector. In 
equation (2b) this adjustment is represented by the 
quantity mD 0 . The analogous correction for the 
derivation in this section, which we will here call d,, 
is given by 

d, = -R K&.r ADJ ) (24) 

which follows from equation (20) by replacing 
(O - Q*) with d, on the left-hand-side and Q with 
D. on the right-hand-side. The D 0 value to be 
inserted into the equations is the most recent value 
determined for operational use. 

Equation (22) can be generalized to allow for 
alignment and/or bias adjustments within the 
calibration state vector. This is done by simply 
combining equations (11a) and (22), with the 
restriction that the summation over k exclude the 
linear scale factor corrections, i.e., 

2£ = -R*K* X B m +R r K* 0 R d 

+ * T X K\ (R K 0 [**' [A 0 R\] D ) [SsJ c (25) 

with the set {Ss k } restricted to nonlinear terms. The 
0 subscript on K 0 and A 0 indicates that the current 
estimate for the gyro alignments is used in construc- 
ting the nonlinear correction coefficients. After the 
calibration set {m, d, {SsJ} has been determined, 
equations (26a) - (26d) can be used to calculate Q. 


n = (G C - D) 

- K . I { (j k , + sa* A 0 (GC- D)] } (26a) 


G = (/, + m) G 0 

(26b) 

D = (/, + m)D 0 + d 


- R.n&.rADj } 

(26c) 

G„ = R„S 1 

(26d) 


Although straightforward to use as the basis for a 
least-squares algorithm (i.e., to solve for the state 
vector {m, d, {& k }} given an error set (CJ). 
equation (25) is somewhat unaesthetic in that it 
mixes a set of parameters pertaining to the com- 
bined gyro system (i.e., {m,d)) with another set 
pertaining to the individual gyroscopes (i.e., {5s k }). 
For elegance in presentation and to support 
engineering analysis of individual gyro behavior, 
having a state vector consisting solely of specific 
parameters of the individual gyros would be 
desirable. Equation (25) could be so recast if we 
were dealing only with sets of three gyros. 
However, for gyro sets containing more than three 
gyros, the parameter set {m, d} captures all of the 
functionally observable information available in the 
maneuver measurements. (Of course, if the full set 
of gyro data is available, the data can be processed 
for each combination of three gyros and the indivi- 
dual gyro parameters extracted, but this defeats the 
processing simplifications discussed herein.) 

This point concerning observability raises a ques- 
tion: for how many gyroscopes can unique scale 
factor information be obtained when equation (22) is 
applied together with the constraint of fixed gyro 
alignments? This question may be readily answered 
for the case where the state vector is restricted to 
linear scale corrections, i.e., In this case, the 8s { 
matrix transforms to an equivalent m matrix via 

m = - K &, A = - [A 1 A] '[A r &, A] (27) 

Both [A T A] and [A T A] can be shown to be 
3-by-3 symmetric matrices, implying that the 
product [A T A] *[A T &, A] is as well. The change 
matrix m therefore has only six independent 
elements, from which we conclude that the 
techniques of this section can provide independent 
scale parameter corrections for at most six gyro- 
scopes. (“At most” applies because any coaligned 
gyroscopes will have degenerate corrections 
irrespective of the total number). 
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Section 5 - Least-Squares Solution and Weight 
Matrix Specification 

For completeness, we present in this section a few 
points pertaining to the selection of weights to be 
applied to the input measurements. As discussed in 
many references on least-squares algorithms (e.g., 
Reference 4), the solution for the batch linear least- 
squares problem associated with a matrix equation 
HX = Y can generally be written as 

X = (H T W H +WJ‘ ( H 1 W Y + W A X A ) (28) 

For our problem, X (the state vector) will be some 
combination of m, d, and/or {Ss t }, Y is [C, T , .... CwT 
for N maneuvers, H is a matrix of state vector multi- 
plying elements constructed from appropriate pieces 
of equation (25), W is a 3N-by-3N weight matrix for 
the error measurements (the elements of Y), and W A 
is a weight matrix associated with the a priori state 
vector estimate, X A . For our problem, because the 
state vector consists of differential changes from the 
previous best estimate, we set X A = 0. Our only 
remaining concern, therefore, is to establish 
reasonable estimates for W and W A . 

Often it is both convenient and reasonable to simply 
set W to 1„ and W A to 0. (We used this approach in 
our analysis of the HST maneuver data and have 
found it operationally acceptable.) Implicit in the 
approach are the following five assumptions: (1) the 
state vector correction terms are fairly stable over 
the time period of data collection, (2) the degree of 
correlation between measurement error components 
is fairly small, (3) the expected error component 
magnitudes are all approximately the same, (4) the 
data set spans the domain of state vector sensitivity 
sufficiently well that observability is not a problem, 
and (5) a sufficiently extensive data set has been 
accumulated that neglect of a priori information 
does not undermine operations. The first three 
points relate to setting W to /, N , whereas the last two 
relate to setting W A to 0. 

If any of the conditions indicated in the previous 
paragraph are significantly violated, a more 
sophisticated weighting scheme is required. We 
present here a method for specifying W that retains 
assumption 1, eliminates assumption 3, and replaces 
assumption 2 with a less restrictive one (called 2a) 
that the measurement errors associated with each 
maneuver are uncorrelated with those of all others. 
We will not consider the possible advantages of a 


nonzero W A . Assumptions 1 and 2a allow W to be 
expressed as a block diagonal matrix, with each 
block being a 3-by-3 matrix, h\ associated with a 
specific maneuver. Given the block diagonal form, 
each w can be written as (p K +p 0 )\ where p„ is the 
covariance associated with reference attitude errors, 
and p a is the covariance associated with random 
gyro errors. The attitude covariance matrix is given 

by 

P, = T; Pi T, + p, (29) 

where p i and p f are the initial and final attitude 
covariance matrices in the instantaneous spacecraft 
frame, and T x is as used in equation (5). Refer- 
ence 5 specifies an equation for attitude covariance 
matrices such as p t and p r This equation, which 
depends upon the reference star distribution and the 
measurement and catalog errors for each star, is 

P« = <*’ [/, - KcXW 7 ] 1 (30) 

where o’ = [Ea/ 2 ] 1 , a, is the root-mean-square 
combined measurement and catalog error for the jth 
star, Vj is the jth star vector expressed in the space- 
craft frame, and the sums are over all observations. 
This expression can be simplified for processing 
purposes in the case of observations from a number 
of well-separated star sensors with fairly narrow 
fields-of-view (narrow relative to the field-of-view 
separations). In this case, each \ i can be replaced 
with the boresight direction vector for the jth sensor 
expressed in spacecraft coordinates, with then 
indicating typical error size for that sensor. This 
substitution eliminates ground processing of the 
reference star data. 

A reasonable, albeit heuristic, model for the 
covariance associated with gyro errors is 

Pa = /,[<vV +<v 1 e 1 ] (3D 

where is the typical single-axis standard 
deviation of the gyro drift rate bias, and is the 
typical scale factor / alignment maneuver error. 
Equation (31) does not attempt to model the physi- 
cal mechanism that produces gyro noise, but rather 
requires the user to provide parameters and 
based on typical spacecraft performance. 
Empirically, for the HST gyroscopes working as a 
set, we find ~ 0.01 arcsecond per second and 
~ 0.2 arcsecond per degree. 
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Section 6 - HST Gyroscope Behavior 

The HST gyroscope system comprises three rate 
gyro assemblies (RGAs) manufactured by 
AlliedSignal Government Electronic Systems. Each 
RGA consists of two single-degree-of-freedom, 
dual-mode, rate integrating, mechanical gyroscopes. 
The high-rate mode has a range of ±1800 degrees 
per hour and a resolution of 7.5 milliarcseconds per 
40-hertz sample; the low-rate mode has a range of 
±20 degrees per hour and a resolution of 
0.125 milliarcsecond per 40-hertz sample. The gyro 
alignments are such that any three can be used to 
completely sample rotations of the spacecraft. The 
onboard system is configured to nominally use four 
gyroscopes simultaneously, keeping the remaining 
two as backups. 

RGA units 2 and 3 (those housing gyros 3, 4, 5, 
and 6) were replaced in December 1993 during the 
first HST servicing mission. All six gyroscopes 
were activated for the servicing mission and early 
on-orbit verification and calibration phase. The 
iterative calibration procedure described in 
References 1 and 6 was followed until convergence 
was achieved. Thereafter, the two gyros in RGA 
unit 1 were deactivated, leaving HST operating with 
four new, freshly calibrated gyroscopes. The active 
gyros are mounted with input-axis unit vectors of 
approximately (±0.586, ±0.617, -0.525), with the 
sign sense for the first two components being 
(--• ++, -+, +-)» f° r gyros 3, 4, 5, and 6, respectively. 
The symmetry of these vectors about the yaw axis is 
significant for understanding the specific mani- 
festation of an observed growing scale error. 

As is typical with spacecraft gyroscopes, the biases 
vary fairly rapidly. For the HST gyroscopes, the 
change in the drift rate bias for both high- and low- 
rate modes has been found to be about 7 arcseconds 
per hour per day. The temporal variation of the 
high-rate mode drift bias vector (i.e., as measured in 
vehicle space) has been found to track the low-rate 
mode vector variations quite closely. This allowed 
implementation of an operational procedure 
whereby only the low-rate mode bias is measured 
frequently, based on data accumulated during 
science pointing with the spacecraft pointing control 
system locked on Fine guidance sensor guide stars. 
The high-rate mode bias is then determined from the 
low-rate mode bias via an additive offset, which is 
monitored for constancy once every 4 to 6 weeks. 
The algorithm used for monitoring the offset had 
been, until recently, essentially that discussed in 
Section 3 in association with equation (13). The 


spacecraft pointing control system is commanded to 
place the gyroscopes in high-rate mode while 
maintaining a constant attitude for approximately 
one orbit (about 95 minutes). Fixed-head star 
tracker star measurements are obtained at the begin- 
ning and end of this constant attitude period and 
used to determine the true attitude change. 

In HST operations, most large maneuvers are pre- 
dominantly about the yaw axis. The predominant 
symptom of the scale factor problem discovered in 
August 1995 was a substantially larger postslew 
pointing error for negative yaw maneuvers than for 
positive yaw maneuvers. Upon examining the 
quantity E = (2£q/0) for maneuvers between the 
time of the first servicing mission and August 1995 
with I t| 3 I > 0.9 and 0 > 90 degrees, we found that 
although the average value of E for positive yaw 
maneuvers stayed near zero, its value for negative 
yaw maneuvers was fairly well fit by the curve 

E - 0.2 + 0.6 ( 1 - e ^ ) 

arcseconds per degree (32a) 

T = 6 months (32b) 

The sense of the error for negative yaw maneuvers 
was such that the spacecraft fell short Of its intended 
destination. The random scatter for E is about 
0.3 arcsecond per degree (3a). 

The analysis techniques described in this paper were 
developed to study the temporal change that was 
seen to have occurred in the HST RGAs. As part of 
our study, we have come to realize that the effects 
of gyroscope nonlinearities are as important as the 
temporal changes that precipitated the study. We 
applied our analysis to a combined set of 
83 maneuvers collected in August 1994 and August 
1995. (Our data indicate that the scale factors had 
stopped changing by August 1994.) For some of our 
analysis runs, we also included a 1-hour period of 
constant attitude. We find that studying the fit 
residuals associated with the constant attitude period 
is important for constructing a high-fidelity model 
of gyroscope response. The results of our analysis 
are specified below. 

(1) To study the change in average linear scale 
relative to the original post-servicing-mission cali- 
bration, we performed a fit using the high mode bias 
offset vector and gyro frame linear scale factors as 
our state vector. The best fit values for this case are 
given in equations (33a) and (33b). 
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doPP.Er = [-1.8X10 1 , 3.4x10’, -7.7xl0 2 f 

± lxlO 2 arcsecond per second (33a) 

[Ss,] c = [5.7x10 s . 4.2x10 s , 

8.4x10 s , 1.74xlOY ± 1x10 s (33b) 

As will be discussed shortly, the bias offset adjust- 
ment is that required to compensate for gyroscope 
nonlinearities, the “true” bias at constant attitude 
already having been eliminated by the standard 
operational procedures. The [5s, ] c elements repre- 
sent the average change in the high-rate mode scale 
factors. The sign sense indicates that the gyros have 
become more sensitive (more counts per degree of 
actual slew). The largest single change, that for 
gyro 6, corresponds to an error of 56 arcseconds for 
a 90-degree slew about the input axis. 

(2) Because of the difference in response for 
positive and negative slews, together with the fact 
that the bias determination procedure had been 
tuned to work accurately at zero angular rate, it 
seemed likely that some scale nonlinearity was 
involved. Taking d = 0 as a constraint effectively 
imposed by the operational procedures, we 
investigated potential nonlinearities by solving for a 
state vector consisting of [5s, ] c and [6s, J c . The best- 
fit results in this case are 

[6s,] c = [6.0x10 s , 2.9x10 s , 

1.27x10"*, 1 .48x10 Y ± lxlO 5 (34a) 

[5s, J c = [0.8x10 s , 6. lxlO 5 , 

1.95X10 4 , 7.8x10 Y ± 1x10 s (34b) 

Comparing the nonlinear correction values with the 
average change values indicated for the first case, 
we see that the error associated with not taking the 
nonlinear effect into account can be as large as the 
temporal change. We also determined fit param- 
eters for two other cases, one including d in the state 
vector and another using g 2 (Q) = £2’ rather than l£2l. 
The former showed a slight reduction in the fit 
residuals, whereas the latter showed a slight increase 
in the fit residuals; the changes in residuals in both 
cases were insignificant. 

Given our findings regarding scale factor non- 
linearities, the spacecraft pointing control logic 
should ideally include compensation for this effect 
when estimating spacecraft angular rates. Although 
the HST pointing control system does not model 
scale factor nonlinearities, we can compensate to a 
significant degree for the nonlinearities by allowing 


the low-to-high bias offset to absorb the average 
effect of the gyroscope nonlinearities as weighted by 
the actual distribution of maneuvers scheduled for 
the HST mission. This is effectively what happens 
with the fit procedure associated with equa- 
tion (33a). The large negative third component for 
the bias in equation (33a) is associated with the 
positive sign of the components of [Ss„] c in equa- 
tion (34b), together with the fact that gyros 3-6 are, 
on average, pointing along the negative yaw axis. 
This weighting for mission maneuver distribution 
will also affect the estimated average scale factors, 
as can be seen by comparing equations (33b) and 
(34a). Empirically, it appears that adequate HST 
mission performance is achieved with this approach 
during normal operations. We note, however, that 
this approach does not give optimized performance 
for high-rate mode, inertial hold conditions, the 
implied spurious drift being about 300 arcseconds 
per hour. 

Using the bias vector to absorb the average effect of 
gyroscope nonlinearities weighted according to the 
profile of mission maneuvers could be problematic 
for spacecraft that use single-mode gyroscopes. For 
such spacecraft, science operations would likely 
require the bias vector to be selected so that pointing 
performance is optimized with respect to constant 
attitude periods. Adjusting the bias to improve 
maneuver performance is therefore not an option. 
Mission engineers designing the pointing control 
and sensor calibration algorithms for such missions 
should consider including compensation for gyro- 
scope nonlinearities, particularly if slewing accura- 
cies better than 1 arcsecond per degree are required. 

(3) As part of our analysis of the HST gyroscope 
changes, we also considered the possibility that the 
changes were associated with the gyroscope 
alignment matrix. We therefore performed a fit for 
a scale factor / alignment correction matrix (m) 
together with a bias adjustment (d) based on 
equation (11a). We found that including the 
alignment adjustments did not significantly improve 
the residuals relative to those associated with the fit 
restricted to state vector [d, [5s,] c }. We specifically 
found that the alignment terms did not allow us to 
simultaneously obtain improved residuals for the 
maneuver data while maintaining small residuals for 
the constant attitude data. Our results are consistent 
with there being no significant change in the 
gyroscope alignments during the 18 months 
following the first HST servicing mission. 
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Conclusions 

This paper has presented a number of variations on 
the Davenport algorithm for gyroscope calibration 
specifically designed to (1) allow analysis with a 
drastically restricted quantity of telemetry data and 
(2) extend the state vector domain to allow study of 
both isolated and nonlinear scale factor corrections. 
We have applied the techniques to data obtained 
during normal operations of HST as part of a study 
of temporal variations of the HST gyroscope scale 
factors. We have found that the HST replacement 
gyroscopes experienced significant change over the 
first 6 to 8 months following the first HST servicing 
mission, the largest individual change corresponding 
to an error in estimated projected rate about the 
input axis of about 56 arcseconds per 90 degrees. 
We have found scale factor nonlinearities that, when 
characterized as differences between scale factors 
associated with positive and negative rotations, are 
as large as 2 parts in 10000, i.e., about 65 arc- 
seconds per 90 degrees. For spacecraft, such as 
HST, that use dual-mode gyroscopes, the effects of 
the nonlinearities can be accommodated to a 
significant degree via adjustments to the high-rate 
mode drift rate bias vector. This approach may be 
inadequate for missions using single-mode gyro- 
scopes. Finally, we find, to within the accuracy of 
our data set, that no significant changes have 
occurred to the gyroscope alignments during the 
first 18 months following the servicing mission. 

The work reported in this article was supported in 
part by National Aeronautics and Space Adminis- 
tration (NASA) contracts NAS 5-31500 and NAS 
5-31000, which enable Computer Sciences 
Corporation and AlliedSignal Technical Services 
Corporation to provide systems engineering, 
analysis, and operations support to NASA’s 
Goddard Space Flight Center. 

Appendix - Model Maneuver Profile 

In this appendix we present the details of one fairly 
common maneuver model. In addition to the total 
maneuver angle, the model uses three input 
parameters characterizing the spacecraft’s maneuver 
execution algorithm. These parameters can be 
selected as the maximum jerk magnitude (J m ), the 
jerk pulse duration (8), and the maximum angular 
velocity magnitude (o^). The maneuver profile is 
symmetric about the midtime (t/2); it is therefore 
sufficient to construct the maneuver profile through 
that time. Throughout the maneuver, the angle (0), 
rate (co), and acceleration (a) are continuous, and the 
jerk (the third time derivative of 0) takes on one of 


three values: J, 0, or -J. The maneuver through its 
midpoint is composed of two, three, or four 
segments, depending upon the value of ©. The con- 
struction for each solution type is presented below. 


Operationally, three auxiliary parameters 
calculated from the three input parameters: 

are first 

- 5 

(A. la) 

0. = 2J„5 3 

(A. lb) 

e b = e.([(e_ 2 +350/28 2 ] +1) 

(A.lc) 


These three equations will be derived below. The 
determination of whether a two-, three-, or four- 
segment half-maneuver pertains depends upon 
where 0 falls relative to 8 t and 0 b ; a two-segment 
solution pertains for 0 in the range [0,8 a ], a three- 
segment solution for the range [0 t ,0J, and a four- 
segment solution for [ 0 b , 7 t]. 

Two-segment solution 

The two-segment solution assumes that the jerk is 
equal to some positive value J for a time period 8 
and equal to -J for a subsequent equal period. The 
functions a(t), co(t), and 0(1) are each required to be 
continuous through the point of discontinuous jerk. 
The angular velocity reaches its maximum value at 
exactly the midpoint of the maneuver, i.e., at 
t/ 2 = 28. The solution for the two segments is 
specified below. 

Segment 1: 0< t < 8 

J(t) = J (J yet unknown) (A.2a) 

a(t) = Jt (A.2b) 

w(t) = 1/2 Jt 2 (A.2c) 

0(0 = 1/6 Jt 3 (A. 2d) 

Segment 2: 8 < t < 28 


J(t) = 

-J 

(A.2e) 

a(t) = 

J 8 - J (i - 8) 

(A.2f) 

co(t) = 

1/2 J 8 2 + J 8 (t - 8) 



- 1/2 J (t - 8) 2 

(A.2g) 

6(t) = 

1/6 J 8 1 + 1/2 J 8 2 (t - 8) 



+ 1/2 J 8 (t - 8) 2 - 1/6 J (t - 8) 3 

(A.2h) 


The unknown J is determined by the requirement 
that 0(t/2) = 8/2. Substituting t = 28 in equation 
(A.2h) yields 

J = e/ 28 3 (A.2i) 
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The two-segment solution applies until equa- 
tion (A.2i) produces a value of J greater than J m . 
This gives the limiting angle 0„ indicated in 
equation (A. lb). 

Three-segment solution 

For maneuvers with angle 0 exceeding ©,, the two 
periods of constant jerk are separated by a period of 
zero jerk, of duration e (to be determined). For 
convenience, let us define a time point A = 8+e. 
The solution for the three segments is specified 
below. 


Segment 1: 0 < t < 8 


J(t) = J n 

(A.3a) 

a(l) = J m t 

(A.3b) 

to(l) = 1/2 J„t 2 

(A.3c) 

0(t) = 1/6 J m t 1 

(A. 3d) 

Segment 2: 8 < t < A 


J(t) = o 

(A.3e) 

a(t) = J m 8 

(A.3f) 

co(t) = 1/2 J m 8 2 + J„ 8 (t - 8) 

(A.3g) 

0(t) = 1/6 J„ 8 3 + 1/2 J m 8 2 (t - 8) 


+ 1/2 J m 8(t - 8) 2 

(A.3h) 

Segment 3: A < t < A+8 


J(t) = -L 

(A.3i) 

a(t) = J m 8 - J m (t - A) 

(A.3j) 

toft) = 1/2 J m 8 2 + J„Se + J m 8(t-A) 


- 1/2 J„(t - A) 2 

(A.3k) 

0(t) = 1/6 J n 8 3 + 1/2 J m 8 2 e + 1/2 J m 

8e 2 


+ 1/2 J m 8 2 (t - A) + J m 8 e (t - A) 

+ 1/2 J m 8 (t - A) 2 - 1/6 J m (t - A) 3 (A.31) 


The unknown e is determined by the requirement 
that 0 (t/ 2) = 0/2. Substituting t = A+8 in equation 
(A.31) yields the quadratic equation 

e 2 + 3 8e - 28 2 (0/0. - 1) = 0 (A.3m) 

the solution for which is 

e = 3/2 8 ([1+8/9 (0/0. - l)] 10 - 1) 

= 1/2 8 {[ l + 80/0J ,n - 3) (A.3n) 

The three-segment solution applies until equa- 
tion (A.3k), combined with equation (A.3n), pro- 
duces a value of to greater than ca^. The maximum 


permitted value of e can be found by setting (o(t) in 
equation (A.3k) to o) m at t = 28+e. This results in 

e n>» = (!)„,/ J„ 8 - 8 (A.3o) 

Note that for the progression of solutions to be 
consistent, we require tq, > J m 8^. The maximum 
maneuver angle permitted for the three-segment 
model can be found by substituting e„„ for e in 
equation (A.3m); the result is equation (A.lc). 

Four-segment solution 

For maneuvers with angle 0 exceeding 0„, the third 
segment is followed by a period of constant angular 
rate at the maximum permitted value. This fourth 
segment lasts until the maneuver reaches the half- 
way point, i.e., until 0(t) = 0/2. The result is that 
the maneuver profile for the first three segments is 
the same as that appropriate for a three-segment 
solution with e = e m , and during the fourth segment 
it is given by 


J(t) = 0 (A.4a) 

a(t) = 0 (A.4b) 

oo(t) = (bu (A.4c) 

0(t) = 0/2 + co m [t - (28+0] (A-4d) 

The total maneuver duration in this case is 
determined by the requirement that 0(t/2) = 0/2. 
Thus, x is given in this case by 

t/2 = (0 - 0J/2O3U + (28+0 (A.4e) 
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